1 Vorbereitung
- 1.1 R-Pakete
- 1.2 Daten: Filmbeurteilung
2 Hingergrund
- 2.1 Bezug zum Studiengang AWM
- 2.2 Forschungsfrage
3 Explorative Analyse
4 Modellefinition
5 Modell in R
- 5.1 Modell 1: Standard-Priors (post1)
- 5.2 Modell 2: Informierte (?) Priors
6 Fazit

1 Vorbereitung

1.1 R-Pakete

library(rstanarm)  # Bayes-Modelle
library(tidyverse)  # Datenjudo
library(bayesplot)  # Plotting
library(gt)  # Tabellen
library(parallel)  # Turbo
library(rstatix)  # Deskriptive Statistiken
library(bayestestR)  # Vernachlässigbare Unterschiede/Zusammenhänge
#library(see)  # Visualisierung
library(tictoc)  # Zeit messen, wie lange das Modell rechnet

Turbo einschalten:

options(mc.cores = parallel::detectCores())

1.2 Daten: Filmbeurteilung

library(ggplot2movies)
data(movies)

Hilfe zu den Daten gibt es hier:

help(movies)

2 Hingergrund

2.1 Bezug zum Studiengang AWM

Nach dem Studium haben Sie bei einem Online-Händler angeheuert und zwar in der Medien-Abteilung. Der Hinweis, dass Sie etwas mit Medien studiert haben genügte. Vielleicht war auch die Behauptung, dass Sie der absolute R-Checker seien nützlich, um den Job zu bekommen … Jedenfalls müssen Sie den Behauptungen Taten folgen lassen, Schluck!

2.2 Forschungsfrage

Wie (stark) ist der Zusammenhang von logarithmierten Budget und Bewertung eines Films?

Unser Kunde – ein reicher Mäzen mit extrentischem Lebenswandel, der gerne ein paar Millionen invetieren möchte – geht von einem positivem Zusammenhang, \(\beta\) aus. Entsprechend sei unsere Hypothese \(H_1: \beta > 0\).

3 Explorative Analyse

3.1 Datensatz vorverarbeiten

Es macht für Analysen, die in Stan laufen immer Sinn

nur die relevanten Variablen in die Analyse aufzunehmen
fehlende Werte zu entfernen
vorab (vor dem Modellieren) etwaige Transformationen (wie Quadrieren, logarithmieren) vorzunehmen

movies <-
  movies %>% 
  mutate(budget_log10 = log10(budget))

movies2 <-
  movies %>% 
  select(budget_log10, rating) %>% 
  drop_na() %>% 
  filter(budget_log10 != -Inf)

Einige Filme haben ein Budget von 0. Das Logarithmus von 0 ist aber minus Unendlich. Mit dieser “Zahl” kann man nicht rechnen. Daher filtern wir alle Zeilen mit -Inf heraus.

3.2 Logarithmus

budget_log10 fasst die Größenordnung des Budgets:

1000: 10^3 -> log10(10^3) = 3
10000: 10^4 -> log10(10^4) = 4
100000: 10^5 -> log10(10^5) = 5

Wir nehmen also die Größenordnung heran, nicht den Betrag selber.

3.3 Deskriptive Statistiken

movies2 %>% 
  get_summary_stats() %>% 
  gt()

variable	n	min	max	median	q1	q3	iqr	mad	mean	sd	se	ci
budget_log10	5183	3	8.301	6.477	5.439	7.176	1.737	1.222	6.232	1.204	0.017	0.033
rating	5183	1	10.000	6.300	5.200	7.200	2.000	1.483	6.137	1.544	0.021	0.042

3.4 Zusammenhang im Datensatz visualisieren

plot0 <- movies %>% 
  ggplot(aes(x = budget_log10, y = rating)) +
  geom_point(alpha = .2)

plot0 +
  geom_smooth(method = "lm")

Es scheint einen sehr schwachen, negativen Zusammenhang zu geben.

Das gefällt unserem Auftraggeber nicht.

Aber, sagt er, und das zu Recht, entscheidend ist ja nicht die Stichprobe, sondern die Population. Weswegen wir uns bitte schön sofort an die Inferenzstatistik bewegen sollten.

Aha, unser Mäzen kennt sich also sogar mit Statistik aus.

4 Modellefinition

4.1 Priors

4.2 Likelihood

\(r_i \sim N(\mu_i, \sigma)\)

mit \(r_i\): Rating für Film \(i\)

4.3 Lineares Modell

\(\mu_i = \alpha + \beta_1 b\)

4.4 Prioris

Die Prioriwerte sind etwas weiter unten aufgeführt.

5 Modell in R

5.1 Modell 1: Standard-Priors (`post1`)

5.1.1 Modellefinition in R (rstanarm)

tic()
post1 <- stan_glm(rating ~ budget_log10,
               data = movies2,
               refresh = 0   # Mit `refresh = 0` bekommt man nicht so viel Ausgabe
               )  
toc()

## 2.042 sec elapsed

Einige Infos zu den Priori-Werten bei stan_glm() findet sich hier.

tic() und toc() messen Anfang und Ende der dazwischen liegenden Zeit. Da die Modelle ein paar Sekunden brauchen, ist es ganz interessant zu wissen, wie lange wir warten mussten.

5.1.2 Priors

Welche Priori-Werte wurden (per Standard) gewählt?

prior_summary(post1)

## Priors for model 'post1' 
## ------
## Intercept (after predictors centered)
##   Specified prior:
##     ~ normal(location = 6.1, scale = 2.5)
##   Adjusted prior:
##     ~ normal(location = 6.1, scale = 3.9)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = 0, scale = 2.5)
##   Adjusted prior:
##     ~ normal(location = 0, scale = 3.2)
## 
## Auxiliary (sigma)
##   Specified prior:
##     ~ exponential(rate = 1)
##   Adjusted prior:
##     ~ exponential(rate = 0.65)
## ------
## See help('prior_summary.stanreg') for more details

Mit coefficients ist das Regressionsgewicht \(\beta\) gemeint.

Das sind keineswegs besonders schlaue Prioris. Könnten wir unseren Mäzen befragen, er scheint ja Experte zu sein, kämen wir vielleicht zu klügeren Prioris (allerdings ist der Mäzen gerade auf einer “Musen-Reise” und nicht zu sprechen).

5.1.3 Posteriori-Verteilung

Überblick über die Parameter:

print(post1)

## stan_glm
##  family:       gaussian [identity]
##  formula:      rating ~ budget_log10
##  observations: 5183
##  predictors:   2
## ------
##              Median MAD_SD
## (Intercept)   6.8    0.1  
## budget_log10 -0.1    0.0  
## 
## Auxiliary parameter(s):
##       Median MAD_SD
## sigma 1.5    0.0   
## 
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanreg

Langfassung:

summary(post1)

## 
## Model Info:
##  function:     stan_glm
##  family:       gaussian [identity]
##  formula:      rating ~ budget_log10
##  algorithm:    sampling
##  sample:       4000 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 5183
##  predictors:   2
## 
## Estimates:
##                mean   sd   10%   50%   90%
## (Intercept)   6.8    0.1  6.6   6.8   6.9 
## budget_log10 -0.1    0.0 -0.1  -0.1  -0.1 
## sigma         1.5    0.0  1.5   1.5   1.6 
## 
## Fit Diagnostics:
##            mean   sd   10%   50%   90%
## mean_PPD 6.1    0.0  6.1   6.1   6.2  
## 
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
## 
## MCMC diagnostics
##               mcse Rhat n_eff
## (Intercept)   0.0  1.0  4026 
## budget_log10  0.0  1.0  4018 
## sigma         0.0  1.0  3182 
## mean_PPD      0.0  1.0  3940 
## log-posterior 0.0  1.0  1765 
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

Nur die mittleren Schätzwerte für die Regression:

coef(post1)

##  (Intercept) budget_log10 
##    6.7632674   -0.1005079

Man kann sich die Posteriori-Intervalle so ausgeben lassen:

posterior_interval(post1) %>% 
  round(2)

##                 5%   95%
## (Intercept)   6.59  6.95
## budget_log10 -0.13 -0.07
## sigma         1.52  1.56

Wir wissen jetzt also schon das Wesentliche: Mit einer Wahrscheinlichkeit von 90% liegt der Zusammenhang (das Beta-Gewicht) für Log-Budget zwischen -0.13 und -0.07. Es gibt also einen gewissen, negativen Zusammenhang zwischen Budget und Bewertung eines Films, laut unserem Golem zumindest.

5.1.4 Visualisieren

5.1.4.1 Priori-Verteilung

post1_prior_pred <- stan_glm(rating ~ budget_log10,
               data = movies2,
               prior_PD = TRUE  # DIESER Schalter gibt uns die Prior-Prädiktiv-Verteilung
               #, refresh = 0
               )  # Mit `refresh = 0` bekommt man nicht so viel Ausgabe

Aus der Hilfeseite:

prior_PD A logical scalar (defaulting to FALSE) indicating whether to draw from the prior predictive distribution instead of conditioning on the outcome.

post1_prior_pred

## stan_glm
##  family:       gaussian [identity]
##  formula:      rating ~ budget_log10
##  observations: 5183
##  predictors:   2
## ------
##              Median MAD_SD
## (Intercept)   5.6   20.4  
## budget_log10  0.1    3.2  
## 
## Auxiliary parameter(s):
##       Median MAD_SD
## sigma 1.2    1.2   
## 
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanreg

Die Koeffizienten aus post1_prior_pred sind also rein durch die Priori-Werte definiert; die Daten sind nicht eingeflossen.

Wenn wir das Objekt mit as_tibble() in eine Tabelle umwandeln, bekommen wir eine Tabelle mit den Stichproben:

post1_prior_pred_draws <- 
  post1_prior_pred %>% 
  as_tibble() %>% 
  rename(a = `(Intercept)`,  # schönere, einfachere Namen
         b = budget_log10) %>% 
  slice_sample(n = 100)

movies2 %>% 
  ggplot() +
  geom_point(aes(x = budget_log10, y = rating)) + 
  geom_abline(data = post1_prior_pred_draws,
aes(intercept = a, slope = b), color = "skyblue", size = 0.2)

Puh, die Priori-Werte sind … vogelwild🐦 .

5.1.4.2 Posteriori-Verteilung: Regressionsgerade

plot1 <- plot0 +
  geom_abline(intercept = coef(post1)[1],
              slope = coef(post1)[2],
              color = "blue")
plot1

col_names <- c("a", "b", "sigma")
draws_m1 <-
  post1 %>% 
  as_tibble() 

names(draws_m1) <- col_names

Ein Blick in die ersten paar Zeilen der Post-Stichproben:

draws_m1 %>% 
  slice_head(n=10) %>% 
  gt() %>% 
  fmt_number(everything(), decimals = 1)

a	b	sigma
6.8	−0.1	1.5
6.6	−0.1	1.5
6.9	−0.1	1.6
6.9	−0.1	1.5
6.7	−0.1	1.5
6.7	−0.1	1.5
6.6	−0.1	1.6
6.7	−0.1	1.5
6.8	−0.1	1.5
6.7	−0.1	1.6

Und hier die Posteriori-Verteilung für \(\alpha\) und \(\beta\) visualisiert:

plot0 +
  geom_abline(data = draws_m1,
              aes(intercept = a,
                  slope = b),
              color = "skyblue1",
              alpha = .1) +
  geom_abline(intercept = coef(post1)[1],
              slope = coef(post1)[2],
              color = "blue")

5.1.4.3 Verteilung von \(\beta\)

Hier die Verteilung für die Steigung (Regressionsgewicht \(\beta\)):

draws_m1 %>% 
  ggplot(aes(x = b)) +
  geom_density()

5.1.4.4 Posterior-Intervalle

Die Posteriori-Intervalle kann man sich schnöde mit plot() ausgeben lassen, wobei man als Parameter den Namen des Modells übergibt.

plot(post1)

Es gibt aber auch andere Darstellungsarten, z.B. als Dichtediagramme:

mcmc_areas(post1) +
  labs(title = "Posteriori-Verteilung",
       caption = "Gezeigt werden Median und 50% bzw. 90% Perzentil-Intervalle")

S. Hilfe hier

Man kann sich auch angeben, welchen Parameter man visualisieren möchte, und zwar mit dem Argument pars (wie paramters).

Oder nur das Posteriori-Interval für den Regressionskoeffizienten:

mcmc_areas(post1,
           pars = "budget_log10")

5.1.5 Fazit

Die Wahrscheinlichkeit, dass der Zusammenhang in Wirklichkeit zwischen -0.15 und -0.05 ist, ist sehr hoch.

Aber was bedeutet das, wie interpretiert man den Befund?

Rufen wir uns dazu ins Gedächtnis, was log10 bedeutet: Es bedeutet, dass man das Budget mit 10 multipliziert, also verzehnfacht, sozusagen eine Null hintendran schreibt.

Also: Verzehnfacht man das Budget, so verringert sich mittlere Bewertung um etwa 0.15 bis 0.05 Punkte.

Ob das viel ist, muss ein Medienexperte beantworten. Vielleicht sind Sie ja einer!

5.2 Modell 2: Informierte (?) Priors

Kramen wir all unser Wissen über Filme und die Filmindustrie zusammen! Wie wichtig sind die Kohlen für die Güte eines Films? Haben die einen großen Einfluss?

(Wenn wir von “Einfluss” sprechen, denken wir sicher automatisch an kausalen Einfluss. Jedenfalls geht es mir so und ich glaube, es ist schwierig, nicht an kausalen Einfluss, sondern nur an blanke Assoization, zu denken).

5.2.1 Modellefinition in R (rstanarm)

Nach langem Beratschschlagen mit Film-Experten gehen wir davon aus, dass es im Prinzip - “normalerweise” keinen Zusammenhang gibt. Allerdings hat man in Filmen schon Pferde kotzen gesehen, deswegen wollen wir unseren Golem zu verstehen geben, dass er sich auch auf Überraschungen einstellen soll …

Kurz gesagt, auf Golem-Sprech:

\[\beta = \mathcal{N}(0,0.2)\]

Mit anderen Worten, verzehnfacht man das Budget, erwarten wir eine Veränderung von (in 95% der Fälle) nicht mehr als ±0.4 Rating-Punkte.

Für die übrigen Priorwerte greifen wir auf wenig ambitionierte Standardwerte zurück.

post2 <- stan_glm(rating ~ budget_log10,
               data = movies2,
               prior_intercept = normal(6, 1),  # alpha
               prior_aux = exponential(1),  # sigma
               prior = normal(0, .2),  # beta
               refresh = 0)  # Nicht so viel Ausgabe

Um Rechenzeit zu sparen, kann man das Modell auch speichern:

save(post1, file = "post2.rda")
load(file = "post2.rda")

5.2.2 Priors

prior_summary(post2)

## Priors for model 'post2' 
## ------
## Intercept (after predictors centered)
##  ~ normal(location = 6, scale = 1)
## 
## Coefficients
##  ~ normal(location = 0, scale = 0.2)
## 
## Auxiliary (sigma)
##  ~ exponential(rate = 1)
## ------
## See help('prior_summary.stanreg') for more details

5.2.3 Prior-Prädiktiv-Verteilung

5.2.3.1 Berechnen der Prior-Prädiktiv-Verteilung

tic()
post2_prior_pred <- stan_glm(rating ~ budget_log10,
                             data = movies2,
                             prior_intercept = normal(6, 1),  # alpha
                             prior_aux = exponential(1),  # sigma
                             prior_PD = TRUE,  # Prior-Prädiktiv-Verteilung
                             prior = normal(0, .2),  # beta
                             refresh = 0)  # Nicht so viel Ausgabe
toc()

## 1.271 sec elapsed

5.2.3.2 Visualisierung der Prior-Prädiktiv-Verteilung

Dazu gehen wir vor wie für Modell 1, post1. Zuerst wandeln wir das Objekt post2 in eine Tabele mit den Stichproben aus der Post-Verteilung um:

post2_prior_pred_draws <- 
  post2_prior_pred %>% 
  as_tibble() %>% 
  rename(a = `(Intercept)`,  # schönere, einfachere Namen
         b = budget_log10) %>% 
  slice_sample(n = 100)

Und dann visualisieren

movies2 %>% 
  ggplot() +
  geom_point(aes(x = budget_log10, y = rating)) + 
  geom_abline(data = post2_prior_pred_draws,
aes(intercept = a, slope = b), color = "skyblue", size = 0.2)

Sieht doch schon viel besser aus. 🏆

5.2.4 Posteriori-Verteilung

5.2.4.1 Kurzfassung

Überblick über die Parameter:

print(post2)

## stan_glm
##  family:       gaussian [identity]
##  formula:      rating ~ budget_log10
##  observations: 5183
##  predictors:   2
## ------
##              Median MAD_SD
## (Intercept)   6.8    0.1  
## budget_log10 -0.1    0.0  
## 
## Auxiliary parameter(s):
##       Median MAD_SD
## sigma 1.5    0.0   
## 
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanreg

Nur die mittleren Schätzwerte für die Regression:

coef(post2)

##  (Intercept) budget_log10 
##   6.75682455  -0.09962434

5.2.4.2 Langfassung

Ausführlicher:

summary(post2)

## 
## Model Info:
##  function:     stan_glm
##  family:       gaussian [identity]
##  formula:      rating ~ budget_log10
##  algorithm:    sampling
##  sample:       4000 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 5183
##  predictors:   2
## 
## Estimates:
##                mean   sd   10%   50%   90%
## (Intercept)   6.8    0.1  6.6   6.8   6.9 
## budget_log10 -0.1    0.0 -0.1  -0.1  -0.1 
## sigma         1.5    0.0  1.5   1.5   1.6 
## 
## Fit Diagnostics:
##            mean   sd   10%   50%   90%
## mean_PPD 6.1    0.0  6.1   6.1   6.2  
## 
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
## 
## MCMC diagnostics
##               mcse Rhat n_eff
## (Intercept)   0.0  1.0  3748 
## budget_log10  0.0  1.0  3789 
## sigma         0.0  1.0  4063 
## mean_PPD      0.0  1.0  3939 
## log-posterior 0.0  1.0  2200 
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

5.2.5 Visualisieren

5.2.5.1 Regressionsgerade

plot1_m2 <- plot0 +
  geom_abline(intercept = coef(post2)[1],
              slope = coef(post2)[2],
              color = "blue")
plot1_m2

Erstellen wir eine Tabelle nur mit den Post-Samples:

col_names <- c("a", "b", "sigma")
draws_m2 <-
  post2 %>% 
  as_tibble() 

names(draws_m2) <- col_names

Die ersten paar Werte aus der Tabelle mit Post-Samples:

a	b	sigma
6.9	−0.1	1.6
6.8	−0.1	1.5
6.9	−0.1	1.5
6.6	−0.1	1.6
6.8	−0.1	1.5
6.9	−0.1	1.5
6.8	−0.1	1.5
6.7	−0.1	1.5
6.7	−0.1	1.5
6.7	−0.1	1.5

Und hier die Regressionsgerade mit dem “Unsicherheitsbereich” für die Regressionskoeffizienten.

plot0 +
  geom_abline(data = draws_m2,
              aes(intercept = a,
                  slope = b),
              color = "skyblue1",
              alpha = .1) +
  geom_abline(intercept = coef(post2)[1],
              slope = coef(post2)[2],
              color = "blue")

5.2.5.2 Verteilung von \(\beta\)

draws_m2 %>% 
  ggplot(aes(x = b)) +
  geom_density()

5.2.5.3 Posterior-Intervalle

S. Hilfe hier

mcmc_areas(post2) +
  labs(title = "Posteriori-Verteilung",
       caption = "Gezeigt werden Median und 50% bzw. 90% Perzentil-Intervalle")

mcmc_intervals(post2,
               pars = "budget_log10")

mcmc_areas(post2,
           pars = "budget_log10") +
  labs(title = "Posteriori-Verteilung",
       caption = "Gezeigt werden Median und 50% bzw. 90% Perzentil-Intervalle")

Das ist das Gleiche wie unser Dichte-Diagramme etwas weiter oben.

5.2.6 Quantile

summary(post2)

## 
## Model Info:
##  function:     stan_glm
##  family:       gaussian [identity]
##  formula:      rating ~ budget_log10
##  algorithm:    sampling
##  sample:       4000 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 5183
##  predictors:   2
## 
## Estimates:
##                mean   sd   10%   50%   90%
## (Intercept)   6.8    0.1  6.6   6.8   6.9 
## budget_log10 -0.1    0.0 -0.1  -0.1  -0.1 
## sigma         1.5    0.0  1.5   1.5   1.6 
## 
## Fit Diagnostics:
##            mean   sd   10%   50%   90%
## mean_PPD 6.1    0.0  6.1   6.1   6.2  
## 
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
## 
## MCMC diagnostics
##               mcse Rhat n_eff
## (Intercept)   0.0  1.0  3748 
## budget_log10  0.0  1.0  3789 
## sigma         0.0  1.0  4063 
## mean_PPD      0.0  1.0  3939 
## log-posterior 0.0  1.0  2200 
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

Laut dem Modell (post2) liegt der Regressionskoeffizient mit 90% Wahrscheinlichkeit eng um -0.1 herum.

Genauer gesagt: \(90\%PI_b: (-0.13, -0.07)\):

draws_m2 %>% 
  summarise(b_90pi = quantile(b, probs = c(0.05, .95)))

## # A tibble: 2 × 1
##    b_90pi
##     <dbl>
## 1 -0.130 
## 2 -0.0702

posterior_interval(post2)

##                      5%         95%
## (Intercept)   6.5688718  6.94726538
## budget_log10 -0.1296204 -0.07015277
## sigma         1.5149537  1.56352142

5.2.7 Wahrscheinlichkeiten für Parameterwerte

5.2.7.1 Positiver Zusammenhang

\(p(b > 0|D)\)

mit “D”, den Daten des Modells.

draws_m2 %>% 
  count(b > 0)

## # A tibble: 1 × 2
##   `b > 0`     n
##   <lgl>   <int>
## 1 FALSE    4000

Die Wahrscheinlichkeit ist praktisch Null, dass der Zusammenhang positiv ist.

5.2.8 \(R^2\)

tic()
post2_r2 <- 
  bayes_R2(post2) %>% 
  as_tibble()
toc()

## 2.2 sec elapsed

post2_r2 %>% 
  ggplot(aes(x=value)) +
  geom_density()

post2_r2 %>% 
  summarise(r2_mean = mean(value),
            r2_median = median(value))

## # A tibble: 1 × 2
##   r2_mean r2_median
##     <dbl>     <dbl>
## 1 0.00623   0.00601

Der Anteil erklärter Varianz ist praktisch Null.

Der Golem blickt es nicht! Das Modell ist schlecht.

5.2.9 PPV

5.2.9.1 PPV berechnen

Simulieren wir den Erfolg neuer Filme; dabei betrachten wir das Budget von \(10^3\) bis \(10^8\) (6 Werte). Wir ziehen pro Budgetwert 1000 Stichproben aus der PPV.

neue_Filme <- tibble(
  budget_log10 = 3:8)

Warum 3 bis 8? Das sind genau die Werte für budget_log10, die wir in daten Daten haben.

Wie viel delogarithmiertem, also “richtigem” Budget entspricht das?

10^(3:8)

## [1] 1e+03 1e+04 1e+05 1e+06 1e+07 1e+08

\(10^3 = 1,000\)
\(10^4 = 10,000\)
…
\(10^8 = 100,000,000\)

ppv_m2 <- 
  posterior_predict(post2, neue_Filme, draws = 1e3) %>% 
  as_tibble() 


dim(ppv_m2)  # Zeilen, Spalten

## [1] 1000    6

Leider werden unsere Eingabewerte (3 bis 8) zerhauen, stattdessen wird 1 bis 6 zurückgegeben. Kümmern wir uns gleich darum.

Hier ein Blick in die Tabelle ppv_m2:

1	2	3	4	5	6
7.176392	6.075872	6.431768	5.281453	5.360870	5.048682
3.460721	6.771258	7.009859	3.730755	5.688797	5.216619
5.637846	5.104807	5.281968	7.426519	4.448275	5.280338
7.500090	8.213477	8.111608	5.958772	7.446054	7.158755
7.382352	6.755573	5.980874	6.033803	7.448571	7.609851
6.442488	6.349062	6.036896	3.925805	3.919113	5.995698

Ändern wir die Spaltennamen von 1,2,…6, in 3,4,…,8 um.

Reine Zahlen werden als Spaltennamen nicht akzeptiert, deswegen wandeln wir nohc die Zahl 3 in den Text "3" um, mit as.character():

names(ppv_m2) <- as.character(neue_Filme$budget_log10)

Vom breiten ins lange Format überführen:

ppv_m2_long <- 
  ppv_m2 %>% 
  pivot_longer(everything(),
               names_to="budget_log10",
               values_to="rating")

Ein paar Erklärungen zu pivot_wider():

ppv_m2_long %>% 
  ggplot(aes(x = budget_log10,
             y = rating)) +
  geom_boxplot()

Tja, keine großen Unterschiede zwischen den Budget-Faktoren (also den Werte von budget_log10). Unser Modell prognostiziert die Daten ganz passabel. Nicht, dass das Ergebnis spektakulär wäre.

5.2.10 90%-Vorhersage-Intervalle

Mit der Funktion predictive_interval kann man sich obige Berechnung sparen, sondern bekommt sie genussfertig nach Hause.

Wir sehen hier die tatsächlichen Rating-Werte pro Budget-Wert, nicht nur \(\mu|b\).

post2_pred <- 
  predictive_interval(post2,
                      newdata=neue_Filme)

post2_pred

##         5%      95%
## 1 3.834045 8.944180
## 2 3.867953 8.868228
## 3 3.676618 8.806842
## 4 3.672329 8.777655
## 5 3.590447 8.533323
## 6 3.443809 8.477426

Wie man sieht, sind die Intervalle sehr groß: Das Modell ist sehr schlecht.

6 Fazit

Die Forschungsfrage war, ob das Budget eines Films mit der Bewertung zusammenhängt.

Dazu wurden zwei einfache lineare Modelle berechnet, die sich in ihren Vorannahmen leicht unterschieden.

6.1 Schätzbereiche für \(\beta\)

Beide Modelle fanden konsistent einen schwachen, negativen linearen Zusammenhang \(\beta\) zwischen Budget und Bewertung: Filme mit mehr Budget wurden konsistent schlechter bewertet, laut den beiden Modellen. Hier sind 90%-PI berichtet:

Modell 1: [-0.13, -0.07]
Modell 2: [-0.13, -0.07]

6.2 Medianer Effekt

Modell 1: -0.1
Modell 2: -0.1

6.3 Beantwortung der Forschungsfrage

Das Modell ist überzeugt, dass es einen leichten, negativen Zusammenhang gibt. Das Modell schließt aus, dass es keinen Zusammenhang gibt.

QM2-Thema5-Fallstudie1: Filmbewertungen

ses

2021-11-09 22:15:55